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Abstract. Analyses of a 3D simulation of the upper layers of a solar convective envelope provide constraints on 
the physical quantities which enter the theoretical formulation of a stochastic excitation model of solar p modes, 
for instance the convective velocities and the turbulent kinetic energy spectrum. These constraints are then used 
to compute the acoustic excitation rate for solar p modes, P. The resulting values are found ~ 5 times larger than 
the values resulting from a computation in which convective velocit ies and e ntrop y fluctuations are obtained with 
a ID solar envelope model built with the time-dependent, nonlocal iGouehl lll977l) extension of the mixing length 
formulation for convection (GMLT). 

This difference is mainly due to the assumed mean anisotropy properties of the velocity field in the excitation 
region. The 3D simulation suggests much larger horizontal velocities compared to vertical ones than in the ID 
GMLT solar model. The values of P obtained with the 3D simulation constraints however are still too small 
compared with the values inferred from solar observations. 

Improvements in the description of the turbulent kinetic energy spectrum and its depth dependence yield further 
increased theoretical values of P which bring them closer to the observations. It is also found that the source 
of excitation arising from the advection of the turbulent fluctuations of entropy by the turbulent movements 
contributes ~ 65 — 75% to the excitation and therefore remains dominant over the Reynolds stress contribution. 
The derived theoretical values of P obtained with the 3D simulation constraints remain smaller by a factor ~ 3 
compared with the solar observations. This shows that the stochastic excitation model still needs to be improved. 

Key words, convection - turbulence - stars:oscillations - Sumoscillations 



1. Introduction 

Solar-type oscillations are believed to be stochastically ex- 
cited by turbulent convection in the near-surface layers of 
the star. The excitation is caused by turbulent convective 
motions which generate acous tic energy which in turn i s 
injected into the p modes (e.g. lGoldreich fc Keelevlll977|) . 
Measurements of the acoustic energy injected into solar- 
like oscillations are among the goals of future space seismic 
missions such as the COROT feaglin fe The Corot Tearr] 
Il998l) and Eddington l|Favata et al.ll2000|) missions. These 
seismic data will make it possible to constrain the the- 
ory of the oscillation excitation and damping, to provide 
valuable information about the properties of stellar con- 
vection, and hence to severely constrain stellar models. 

Models for stochastic excitation of stellar p modes 
have been proposed by several authors, (e.g. 
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semianalytical approaches yield the acoustic energy 
injected into solar-like oscillations. This offers the advan- 
tage of testing separately several properties entering the 
excitation mechanism which are not well understood or 
modeled. 

Such approaches require simplifying assumptions 
which need to be validated before they can be used with 
confidence. They require an accurate knowledge of the 
properties of turbulent convection and, unfortunately, cur- 
rent observations of the solar granulation cannot provide a 
determination of the turbulent spectrum precise enough in 
the p resent context l|Rieutord et aljEuOOHNordlund et alJ 
1997). On the theoretical side, theoretical models of turbu- 
lent convection, such as the mixing- length approaches or 
multiple size eddies appro aches fe.g. lCanuto fc MazzitellJ 
199ll ICanuto et aHll9 96h provide a too limited descrip- 
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tion of the characteristic scale length of the solar turbulent 
spectrum. 

These theoretical formulations of stochastic excitation 
also involve scaling parameters which are determined to 
by the requirement that the computed values of the oscil- 
lation amplitude s give the best fit to the solar seismic mea- 
surements (e.g. iHoudek et all Il999t ISamadi et all 120011 
Paper II hereafter). When the scaling parameters are so 
adjusted, constraints and validation on the turbulent stel- 
lar medium can only come from seismic observations of 
other stars. Such accurate data on the excitation rates for 
other stars than the Sun are not yet available. 

An alternative way is then to consider results from 3D 
numerical simulations. They indeed enable one to compute 
directly the rate at which p modes are excited (e.g. th is 
was undertaken for the Sun bv lStein fc Nordlundl l)200ll) 1. 
Such methods are time consuming and do not easily allow 
massive computations of the excitation rate for stars with 
different temperatures and luminosities. They can provide 
quantities which can be implemented in a formulation for 
the excitation rate P. In any case we cannot avoid to use 
a ID model for computing accurate eigenfrequencies for 
the whole observed frequency range. 

The purpose of the present paper is to provide a bet- 
ter insight into the excitation model with a semianalytical 
approach but using a model of turbulence and values of 
the scaling parameters derived from a 3D simulation of the 
solar outer layers. We consider in this w ork the theoretical 
formu lation of stochastic excitati on byj Samadi & Goupil 
(2001, hereafter Paper I, see also lSamadil l)200l|) for a de^ 
tailed summary) which includes a detailed treatment of 
turbulent convection. This formulation involves two scal- 
ing parameters which are related to the spatial and tempo- 
ral characteristics of the turbulence model. Our final goal 
is to test the excitation model without adjusting these 
parameters and without the use of the mixing-length ap- 
proach for estimating convective velocities and entropy 
fluctuations. 

The paper is organized as follows: in Sect [21 we briefly 
recall the adopted formulation for estimating the rate at 
which turbulent convection supplies energy to the p modes 
(excitation rate P{v)). We emphasize some assumptions 
and approximations entering this formulation. 

In Section |21 a 3D numerical simulation of the upper 
part of the solar convection zone is used in order to de- 
termine the time averaged properties of turbulent convec- 
tion: this provides constraints on the ingredients involved 
in the theoretical expression of the excitation rate, such 
as scaling parameters, velocity anisotropy factor, the val- 
ues of convective velocities and entropy fluctuations and 
the k (wavenumber) dependence of the kinetic turbulent 
spectrum. 

These constraints are then used in Section 01 to com- 
pute the excitation rate P(v), for radial solar p modes. 
The results are compared with solar seismic observa- 
tions as given in IChaplin et ID ilQOflh and wi th a ID 
mixing-length model built according to iGougbl l)l977l) 's 



non-local formulation of the mixing-length theory (GMLT 
hereafter). In Section |5l we summarize our results and 
discuss some possible origins of the remaining discrep- 
ancies with solar seism ic observations and results by 
IStein fc NordJlElOll). 

2. Stochastic excitation 

2.1. The excitation model 

The rate at which turbulent motions of the convective 
elements supply energy to acoustic oscillation modes is 
computed as in Paper I. For a given mode with eigenfre- 
quency luq, the excitation rate can be written as (Eq. 58 
and Eq. 59. of Paper I) : 



P(w ) =Pr + Ps 



where 



M 



Pr,s = 2j / dmp a — wF R 



(1) 



(2) 



where po is the mean averaged density, w is the rms value 
of the vertical component of the velocity, 



I'M 

1=1 dmg 
Jo 



(3) 



is the mode inertia, £ r is the radial component of the 
fluid displacement adiabatic eigenfunction and $ is an 
anisotropy factor. Following iGoughl l)l977ft . we define 



*(*) 



< M z > - < U > z 



(4) 



where u is the velocity field, < . > denotes horizontal 
average and () denotes time average. The mean vertical 
velocity, w, is defined as: 



w 2 = < ul > - < u, > 2 



(5) 



Pr-, Ps respectively account for the excitation by the 
Reynolds stress and for the excitation resulting from the 
advection of the entropy fluctuations by the turbulent ve- 
locity field (the so-called entropy source term). Here the 
entropy term (Fg) is an advective term which mixes tur- 
bulent pressure and entropy fluctuations. Expressions for 
F R , F s are : 



Fr = f R (0 S r (lj ) F s = f s (0 S s {lj ) 
with 



fs(0 = X 



a s s \ 2 g r (£ r ,ra) 
p Q w 



(6) 

(7) 
(8) 



where a s = (dp/ds) , p denotes the pressure and s the en- 
tropy, s is the rms value of the entropy fluctuations, and 
S and "K are anisotropy factors. We assume that injec- 
tion of acoustic energy into the modes is isotropic. This 
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assumption implies S = 16/15 and % = 4/3 in Eq.Q and 
Eq.JSJ) above. Effects of the space averaged anisotropy in 
the driving process is investigated in Sect. I4.1T21 
The function g r (£ r , m) is denned as : 



stochastic excitation occurs. The gaussian function takes 
the form 



1 



-{ui/ui k ) 



(13) 



, 1 da s d£ r 

g T {£ T ,m) = — — 

a, ar ar 



dr 2 



(9) 



One can show that the Reynolds contribution - Pr - 
scales as $ 2 u> 4 while the source term involving the entropy 
fluctuations - Ps - scales as $ w 2 s 2 . P(uj) is thus very de- 
pendent of the estimated values of w 2 , s 2 and $. The MLT 
provides estimates for w but $ i s a free parameter. For 
isotropic turbulence $ = 3, and in lBohm - Vitensel 1 1958L 
BV-MLT hereafter) formulation <I> = 2. In the present pa- 
per, unless otherwise stated, $ is given by a simulation 
of the upper part of the solar convective zone in Sect. 13.41 
below. 



where Lo k is its linewidth. 

Let Tk be the characteristic time correlation length of 
an eddy of wavenumber fc. Eq. (|13|) corresponds in the 
time domain to a gaussian function with linewidth equal 
to 2/wfc. Then uj k = 2/r k for a gaussian time spectrum. 

The energy supply rate P crucially depends on the cor- 
relation time-scale r k (see Paper II). Following iBalmforthl 
l|l992a[) we define it as : 



For the driving sources in Eq. © 



St 



dk E(k,r) 



E(k,r) 



u re u 

dk E(k,r) E s (k,r) 
k 2 Uq s 2 

\ / dw»(wo+w)xt(w) 



u, 



II 



Ss = 



(10) 



(11) 



E{k) represents the kinetic energy spectrum associ- 
ated with the turbulent velocity field and E s {k) models 
the spectrum of the turbulent entropy fluctuations, with 
k the eddy wavenumber. The time-dependent part of the 
turbulent spectrum is described by the function Xfc(^) 
which models the correlation time-scale of an eddy with 
wavenumber k. The quantity uq = y$/3w is introduced 
for convenience (see Eq. I17J) , 

The above expression for P is mainly based on the 
assumption that the medium is incompressible. In other 
words, we adopt the Boussinesq approximation i.e. assume 
a homogeneous model for the turbulence and the excita- 
tion mechanism. We therefore neglect effects of the strat- 
ification in the excitation process. 

2.2. The turbulence model 

Let ko(r) be the wavenumber at which energy is injected 
into the turbulent cascade and the energy E(k) is maxi- 
mum, fco (r) is related to the mixing- length A = a H p by 
(Paper I): 

2tt 2tt 



„MLT 



(r) 



(12) 



f3A(r) f3aH p {r) 

where (3 is a parameter of order unity, a is the mixing- 
length parameter and H p is the pressure scale-height. This 
is a natural way to estimate kg as A is the characteristic 
length of the largest convective elements. 

The ga ussian function is usually assumed for modeling 
Xfe(w) (e.g. lSteinll967tlGoldreich fc Keelevll977j) as a con- 
sequence of the turbulent nature of the medium where the 



r k = A (ku k y 



(14) 



where u k is the velocity of an eddy with wavenumber k. 
The velocit y u k is obta ined from the kinetic energy spec- 
trum E(k) l|SteirJll967|) 



2/. 



dkE(k). 



E(k) is normalised such that 

f°° 1 

/ dfc E(k) = -<u 2 >-<u> 
Jo 2 



7 = ^<{z) 



(15) 



(16) 



where Uq is introduced for convenience. According to 
Eq. (0} and (|16fl . uq and w are then related to each other 

by 



\u 2 Q = ^{z)w 2 {z) 



(17) 



The parameter A in Eq. (|14fl accounts for our lack of 
precise knowledge of the time correlation r k in stellar con- 
ditions. In the present paper, we assume A = 1 while (3 
fEq. H2fl and $(z) (Eq.^J are given by a simulation of the 
upper part of the solar convective zone in Sect. 13.41 below. 

2.3. Computations of the excitation rate P{uo) 

In practice, we compute the excitation rate P(lu) accord- 
ing to Eq. The calculation requires the knowledge of 
several quantities which can be obtained either from a 
ID model (Paper II) or at least partly from a 3D simula- 
tion. Comparison of the results using both options yields 
insights in the excitation mechanism and its modelling. 
Hence in the following: 

• The velocity, entropy fluctuations, anisotropy and 
turbulent spectra E' s are obtained from a 3D simulation 
as described in the next section. 

• The mean density, po, the thermodynamic quantity 
a s , the oscillation properties - eigenfrequencies and eigen- 
functions - are calculated from a solar envelope equilib- 
rium model and IBalmforthl l)l99 2b Ts pulsation code. The 
envelope model is built with a treatment of convection 
as prescribed by t he GMLT formulati on an d is computed 
in the manner of IBalmforthl l|l992bj) and iHoudek et all 
( 1999J). This solar envelope model (hereafter GMLT solar 
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mode l) is identical to the one considered in l|Samadi et alJ 
120021 hereafter). In particular, it incorporates turbulent 
pressure (momentum flux) in the equilibrium model enve- 
lope. The entire envelope is integrated using the equations 
appropriate to the nonlocal mixing-length formulation by 
Gough (1976) an d to the Eddington ap proximation to ra- 
diative transfer ljUnno fe Spiegell 179661. The equation of 
state included a detailed treatment of the ionization of C, 
N, and O, and a treatment of th e ionization of the next 
seven most abundant elements l)Christensen-Dalsgaardl 
as well as 'pressure ioniz ation' by the metho d of 
Eggleton, Faulkner & Flannery l|Eggleton et alJll973[) . In 
this generalization of the mixing- length approach, two ad- 
ditional parameters, namely a and b, are introduced which 
control the spatial coherence of the ensemble of eddies con- 
tributing to the total heat and momentum fluxes (a) , and 
the degree to which the turbulent fluxes are coupled to the 
local stratification (b). These convection parameters are 
calibrated to a solar model to obtain the helioseismically 
inferred depth of the solar convection zone of 0.287 of the 
solar radius (Christensen-Dalsgaard, Gough & Thompson 
1991). The adopted value for the shape factor $ = 1.3745, 
a value which provides the best fit be tween computed solar 
dam ping rates and measu rements bv lChaplin et alJ l)l998|) 
(see iHoudek et aLlEoOlh . The detailed equations describ- 
in g the equilibrium a nd puls ation models w ere discussed 
bv lBalmforthI <ll992bh and bv lHoudekl lll99^ . 



3.1. Fourier transforms and averaging 

We compute the 2D Fourier transform, along horizontal 
planes, of the velocity field u and the entropy s, at each 
layer z. This provides u(k, z, t) and s(fc, z, t) where k is the 
wavenumber along the horizontal plane. Next we integrate 
u 2 (k, z, t) and s 2 (k, z, t) over circles with radius fc at each 
given layer z. Finally take a time average of the various 
quantities over the time series. This yields u(k, z) and 
s(k, z) where k = \\k\\ is the wavenumber norm. 

We define the time averaged kinetic energy spectrum 
E(k, z) as: 



E(k,z) 



\ u 2 (k, z) for k > 



(18) 



for fc = 

and the time averaged spectrum of the entropy E s (k,z) 



E s (k lZ ) = \ 2 



\s 2 {k,z) for fc > 



for fc = 



(19) 



From Parseval-Plancherel's relation, E(k, z) and 
E s (k, z) satisfy: 



+00 



dkE(k,z) 



<u* > - <u> 



dkE s (k,z) 



2 = 2 UI{Z) 
- < s 2 > - < s > 2 = - s 2 (z) 



(20) 



For implementation in Eas. (|llll(l . the quantities from 
the 3D simulation are interpolated at the GMLT model 
mesh points. The grid of mesh points of the simulated 
domain is matched with the GMLT one such that w in 
the 3D simulation has its maximum at the same layer as 
in the GMLT model. In the simulation, w peaks ~ 40 km 
above the layer at which the mean optical depth < t > is 
unity while in the GMLT model, w peaks ~ 130 km below 
the photosphere (< r >= 2/3). 



3. Constraints from the 3D simulation 

We consider a 3D simulation of the upper part of the solar 
convective zone obtained with the 3D numerical code de- 
veloped at the Niels Bohr Institute for Astronomy, Physics 
and Geophysics (Copenhagen, Denmark). 

The simulated domain is 3.2 Mm deep and its surface 
is 6 x 6 Mm 2 . The grid of mesh points is 256 x 256 x 
163, the total duration 27 mn and the sampling time 30s. 
Physic al assumptions are described in lStein fc Nordlundl 
(1993). 



Output of the simulation considered here are the ve- 
locity field u(x, y,z,t) and the entropy s(x,y, z,t). They 
are used to determine the quantities s 2 , < &(z), w(z), 
E s (k,z), E(k,z) which enter the excitation rate through 

Eqs. BEBHDH5. 



Hence the definitions of the energy spectra here involve 
zero mean velocity and entropy fluctuations. 

3.2. Convective velocities and entropy fluctuations 

Figs, n 121 present w(z), s 2 (z) versus depth for the 3D 
simulation. For comparison purpose, the plots also show 
w and s 2 obtained with the GMLT solar model. 

The vertical velocity GMLT w is larger at the top of 
superadiabatic region but smaller just beneath compared 
to values from the simulation. The GMLT s 2 is larger 
than in the simulation (~ 20 %). This explains that the 
relative contribution of the entropy source term to the 
excitation is overestimated with the GMLT model (see 
Sect. ~ 



14.1.10 . Differences between the GMLT and the sim- 
ulation are likely to be related to differences in the convec- 
tive efficiency: G MLT is less efficient than the 3D. Indeed 
as pointed out bv IHoudek fc Goughl ll2002h . a single eddy 
approach such as the GMLT results in a larger peak for 
the superadiabatic gradient. 

3.3. Velocity anisotropy at large scale 



As it will be shown in Sect. 14.1.21 the value of $ plays 
a crucial role in controlling the depth of the excitation 
region and therefore the total amount of acoustic energy 
injected into the oscillation modes. 

Fig. displays the anisotropy factor <!> versus depth 
z for the 3D simulation. <J>(z) sharply decreases from the 
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Fig. 1. The root mean square of the vertical component 
of the velocity (w = y/ < u 2 . > — < u z > 2 ) in the upper 
layers of a solar model is plotted versus depth for the 
3D simulation (solid line) and for the ID GMLT model 
(dashed line). The abscissa is the depth z = r — i? Q where 
i?0 is the radius at the photosphere. The w maximum 
corresponds to the top of the superadiabatic region and is 
reached at the depth z ~ —130 km in the GMLT model. 
The grid of mesh points of the simulated domain is ad- 
justed so that the w maxima of the 3D simulation and the 
GMLT coincide at the same layer. 




- R„ [Mm] 



Fig. 2. Same as Fig.^for the mean square of the entropy 
fluctuations (s 2 ). The peak is narrower than for the veloc- 
ity w because s 2 scales approximatively as w A . 



value $ = 3 at the top of the CZ down to <E> = 2 and then 
slowly decreases to reach the value $ ~ 1.3 at the bottom 
of the simulation. The decrease of $(z) with depth is ex- 
plained first by the onset of the convection and the forma- 
tion of convective plumes at z ~ and then by the relative 
increase in number of the plumes inward in the simulation. 
Indeed, plumes are highly anisotropic structures whereas 
turbulent cells are quite isotropic. The turbulent Mach 
number increases with z and reaches its maximum value 
at the top of the CZ. Therefore the fluid is more turbu- 



Fig. 3. Same as Fig. ^ for the anisotropy factor <E> versus 



lent outward in the atmosphere. Consequently the number 
of turbulent isotropic cells increases with z up to the top 
of the CZ whereas the number of plumes remains roughly 
constant. The medium is thus more isotropic outward than 
inward. 

In most of the excitation region, the value of $ = 2 
consistent with the BV-MLT is in better agreement with 
the values of $(z) inferred from the simulation compared 
to the value $ = 1.3745 which must be imposed for the 
GMLT solar model in order to match the observed solar 
damping rates. 

3.4. Turbulent kinetic energy spectrum E(k) 

Variations of E and E s with fc at different depths z are 
depicted in Fig. The spectra clearly show two regimes: 
at large scale (small values of fc), the spectra increase ap- 
proximately as k +1 which can roughly be explained using 
dimensional analysis. At small scale (large values of fc), the 
spectra decrease very rapidly with fc. The Kolmogorov law 
(fc -5 / 3 ) is observed only over a small fc-range. Departures 
of the computed spectra from a Kolmogorov law at high 
values of k can be explained by the finite resolution of the 
simulation spatial grid. 

The main characteristics of the kinetic spectrum 
E(k, z) - k dependency - derived from the 3D simula- 
tion are approximatively reproduc ed by an analyt i cal ex - 
pression which was considered bv iMusielak et all l)l994|) . 
namely the 'Ext ended Kolmogorov Spe ctrum' (EKS here- 
after) defined in lMusielak et al.1 1 1994h as: 



E(k,z)=a-§ { 



+ 1 



for fc < fc<f (z) 



-5/3 



(21) 



for fc > fc(f (z) 
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where a is a normalisation factor which satisfies Eq. (|20[l 
and uq is defined according to Eq. l(2Ujl . fc(f is the scale of 
maximum energy in the energy spectrum. 

At each layer, is determined by imposing that the 
EKS, as defined above, matches the turbulent spectrum 
E(k, z) calculated from the simulation as well as possible. 
This then fixes the z dependency of k^ . A similar proce- 
dure is applied for E s {k, z) for which we introduce k^ s . All 
spectra satisfy their respective normalisation condition as 
given in Eq. (J20J . 

For comparison, in Fig. [7| the 'Nesis Kolmogorov 
Spectrum" ( NKS hereaf t er) d etermined from solar ob- 
servations of iNesis et alJ ill 993j) is also shown. The NKS 
scales as fc -5 in the energy injection region for k < k^ 
and down to k min = 0.7 . This spectrum does not agree 
with turbulent spectrum E(k, z) calculated from the sim- 
ulation. In particular, the NKS underestimates the veloc- 
ity of the small size turbulent elements in the cascade 
(k > fco) and overestimates the velocity of the turbulent 
with wavenumber fc ~ fcjf . As we will show in Sect. 14.21 
differences between the EKS and the NKS have an impor- 
tant impact on P(ui). 

If we assume that fejf = k® s , one can show that P(uj) 
scales as fcg~ 4 . P(oj) is therefore very dependent on the 
values reached by fco (z) in the excitation region. Variation 
of fco (2) with depth is thus shown in Fig. for E and E s ; 
fcg 5 " and fcg 5 vary slowly within the excitation region. 

For comparison, in Fig.[S]we have also plotted fc VILT (2;), 
the MLT value for fco (z) according to Eq. (|12|l . The scaling 
parameter j3 in the definition of fc VILT is determined such 
that k^ 111 and fcp take the same value at the layer z ~ 
— 130 km where w reaches its maximum (and consequently 
the layer where the excitation is maximum). The derived 
value is /3 = 3.48. 

fc^LT varies slowly with depth below the top of the 
superadiabatic region (z ~ —130 km) but increases very 
rapidly above. Such a behavior is explained by the rapid 
decrease with z of the pressure scale height H p (which en- 
ters in the definition of fc^ 1 ", Eq. (|12fl ) in the atmosphere. 

Comparison between fc VILT (z) and fcjfO 2 ) shows that 
the mixing-length approach does not model satisfactorily 
the behavior of k${z) in particular just above the layer 
at which w reaches its maximum value. Consequences in 
terms of mode excitation are investigated in Sect. FOl 



4. Consequences in term of p modes excitation 

The acoustic energy supply rate P injected into the solar 
oscillations is related to the rms value v s of surface velocity 
as : 



P(w ) = 2 ?/ 



(22) 



where 77 is the mode damping rate and r s is the radius at 
which oscillations are measured. 

We derive the 'observed' P from IChaplin et alJ 
l)l998[) 's seismic data according to Eq. i\T3 where the 
mode damping rate ,77, and the mode surface velocity , 
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Fig. 4. Turbulent kinetic energy spectra E (top) and E 8 
(bottom) from the simulation are plotted versus fc and 
for different depths z in the simulation. The straight solid 
lines delimitate the slopes fc 1 and fc -5 / 3 of the NKS spec- 
trum. Intersection of the slopes determines fcjf, the scale 
of maximum energy at each depth z 



7j s , are obtained from lOianlin et alJ il99Sft 'a data. The 
mode mass I/^(r s ) is given by the GMLT model and we 
adopt r s = Rq + 200 km consistent with the observations. 

Theoretical values of P are computed according to 
Eq. QJ. In Eas. (|10Jl and (|11|1 the integrations over fc are 
performed from fc = fc m i„ (where fc m m depends on the 
adopted turbulent spectra E and E s ) to fc = 20 fco- We 
checked numerically that contributions to the excitation 
rate from turbulent elements with fc > 20 fco are negligible. 

A gaussian function is assumed for Xfc( w ) m Eq. lflT)|l 
and Eq. fTTj) . 

For the other quantities (w, s 2 , 4>, fco, E(k/ko) and 
E s (k/ko)) involved in the expression for P we investigate 
several possible assumptions. 



4.1. Convective velocities and large scale anisotropy 

In this section, the excitation rate P (Eq. (JJJ) is computed 
with the following assumptions: 
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Fig. 5. The wavenumbers k$ (solid line) and k " (dashed 
line) are plotted versus z (fejf and fcjf 8 are obtained by fit- 
ting, at each layer, the EKS fEq. ETj l to the computed spec- 
trum E and E s of Fig. 0]resp., see text for details). The 
dotted line corresponds to k^f LT (z) obtained according to 
Eq. JX2). In computing k^ LT (z), we assume (3 = 3.48 in 
order for fc^ LT to match the value reached by (solid 
line) at the layer z ~ —130 km where w reaches its maxi- 
mum (/e,f = 3.62 Mm -1 at that layer). 



- the ^-dependency E and E s is given by the analytical 
form of Eq. (2TJ, also called the EKS. 

- k$ = fc(f s = k^ LT where k™ LT (z) is given in Eq. 
with — 3.48 so that fc^ LT takes the value reached by 
kg(z) (solid line, Fig.[5J) at the layer z ~ —130 km where 
w reaches its maximum. 

For the quantities <&, w and s 2 we investigate the effects 
of using either the values derived from the 3D simulation 
(see Sect. El and Sect. EH) or calculated with the GMLT 
solar model. 

4.1.1. Convective velocities and entropy fluctuations 

The values of w, s 2 and $(z) are fixed by the 3D simu- 
lation inside the simulation domain and by the ID equi- 
librium model outside this domain. Either, if we impose 
zero values or if we assume quantities from the ID MLT 
model, no sensitivity on the calculation of P is found. 

Results are shown in Fig. [fj] for P and for the relative 
contribution of the Reynolds stress to the total energy 
supply rate P. When the excitation rate P{v) is com- 
puted with quantities derived from the 3D simulation as 
described in Sect. 2. 3, the resulting excitation rate at max- 
imum is found too small by a factor ~ 4 compared with 
the observations. 

Provided the appropriate value for $ is given in the 
GMLT estimations (see Sect. 4. 1.2 below), no significant 
difference is found in the excitation rate when computed 
with the values of w and s 2 from the simulation or their 
respective GMLT estimations. 



The main effect is illustrated in the bottom panel of 
Fig. EI the 3D simulation generates a larger relative contri- 
bution of the Reynolds stress to P than the GMLT model. 
This is explained as follows: within most part of the exci- 
tation region - except at the top of superadiabatic region - 
the values reached by w are larger whereas values reached 
by s 2 are smaller than their corresponding GMLT estima- 
tions. 

4.1.2. Velocity anisotropy at large scales 

The main consequence (in term of p modes excitation) of 
the differences between the time averaged properties of 
the convective region inferred from the 3D simulation and 
from the GMLT solar model (Fig. |SJ is due to differences 
in their respective anisotropy factor $ values (Fig.[3J). 

Within most of the excitation region, $(z) is found 
close to ~ 2 and thus larger than the value $ = 1.37 as- 
sumed for the ID equilibrium model (see Sect. 13.31 and 
Fig. El . Smaller values of <& decrease the rms total con- 
vective velocity which results in larger values of (see 
Eq. 1141 Eq. I15J) and therefore in a smaller depth of exci- 
tation for a given mode frequency (see Paper III for more 
details). Smaller values of the rms total convective ve- 
locity also induce smaller values of E in the integrand of 
Eq. Q). Consequently, as it is illustrated in Fig.[f)| the total 
amount of acoustic energy injected into the modes is ~ 5 
times smaller for the constant value $ = 1.37 compared 
to the constant value 4> = 2 (the relative contribution of 
the Reynolds stress to P is found ~ 2 times larger in the 
simulation) . 

The effect of the depth dependency of $ on the mode 
excitation is small except at high frequency. This is illus- 
trated in the bottom panel of Fig. H3 (compare the solid 
line with the dashed line). Just above the top of the su- 
peradiabatic region (z > —130 km), $>(z) increases rapidly 
with z until the value ~ 3 (Fig. 0J. Most of the injection 
of acoustic energy into the high frequency modes occurs 
at the top of superadiabatic region. The high frequency 
modes are therefore more sensitive to this rapid increase 
of &(z). As a consequence, the relative contribution of 
the Reynolds stress is larger for the high frequency modes 
than it is when assuming the constant value $ = 2. 

4.2. Turbulent spectra 

In this section we compare the excitation rate obtained 
assuming, for the turbulent spectra (E and E s ) either 
- the EKS spectrum (Eq. |2J with slopes given by the 
3D simulation as in the previous section- or assuming 
the N KS spectrum from solar observations of iNesis et alJ 
l|l993|) (see also Paper I). 

As in Sect. l4.ll we compute P using w, s 2 and <& derived 
from the 3D simulation and assuming that k$ = k^ s = 
£^ LT . The results are plotted in Fig. El 

The NKS overestimates the maximum in P by a factor 
~ 1.5 while the EKS underestimates it by a factor ~ 4. 
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Fig. 6. Top: Rate P at which acoustic energy is injected 
into the solar r adial modes. The fil led dots represent P 
computed from IChaplin et all (|l998f> 's solar seismic data 
according to Eq. Ij22(l . The curves represent theoretical 
values of P computed according Eq. and for differ- 
ent computations of w, s 2 and <3>(z): Solid line: values of 
w, s 2 and $(z) are fixed by the 3D simulation inside the 
simulated domain and by the ID equilibrium model out- 
side this domain. Dashed line (— — — ): same as solid 
line but fixing $ to BV's value <& = 2. Dot dashed line 
(— . — .— ) : values of w, s 2 and $ (= 1.37) are fixed 
by the ID equilibrium model (GMLT). Three dots dashed 
line (— ... — ...—): same as the dot dashed line but 
fixing $ to BV's value $ = 2. Bottom: Same as top 
panel for the relative contribution of the Reynolds stress, 
Pr, to the total acoustic energy P. 



This is because most of the kinetic energy in the NKS is 
concentrated at k ~ k^ whereas in the EKS a large part 
of the energy is concentrated both at large scales (k < k^ ) 
and at small scales (k > k$). 
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Fig. 7. The NKS and the EKS turbulent kinetic energy 
spectra are plotted versus the normalized wavenumber 
k/ko- 
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Fig. 8. Acoustic energy supply rate P computed accord- 
ing to Eq. and assuming for E(k) the EKS and the 
NKS plotted in Fig. Dots represent the energy supply 
rate injected into the oscillations derived from the solar 
observations with the help of Eq. I|22|l . 



4.3. Effects due to the stratification of the turbulent 
spectrum at large scale 

In Sect. 13.41 we showed that the variations of k$ and k^ 3 
with z deduced from the 3D simulation differ from the 
MIT estimation as given by Eq. JT2J| (see Fig. [5}. Fig.|§l 

presents the consequences of the z variations of on 
the oscillation amplitudes (as the variations of fcjf and 
kg 3 with depth are quite similar we assume for the sake 
of simplicity that fc(f s (z) is equal to k^iz)). The z de- 
pendency of fcg causes the maximum of P to be larger 
than when assuming = fc^ LT (~ 50 % larger) . This is 
due to the fact that in most part of the excitation region 
k^ 11 ^ is smaller than k$ and k^ 3 except above the top of 
the superadiabatic region (see Fig. |SJ) . A larger fcjf results 
in a larger linewidth ujk = (kuk)/\ for Xfe( w ) hence in a 
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larger amount of acoustic energy injected to the mode (see 
Eq.EHJl. 

Furthermore, at high frequency P decreases with v 
more rapidly than when assuming k® — k^ s — kJf~ LT . 
Taking into account the actual variation k$ with z instead 
of assuming k^ — fc^ LT makes then the ^-dependency of 
P at high frequency closer to that of the observed exci- 
tation spectrum. This is because, above the top of the 
superadiabatic region, k$ decreases with z whereas fc^ LT 
increases with z . Indeed, the excitation of the high fre- 
quency modes occurs predominantly in the upper most 
part of the top of the superadiabatic region. As men- 
tionned above the line width of Xfct^) decreases with de- 
creasing kg. Therefore the contribution of the term Xfe(w) 
to the excitation of high frequency mode is smaller when 
assuming the actual variation of kg with z than when 
assuming that k® varies as fc^ LT . 



5. Summary and discussion 

An analysis of a 3D simulation of the upper part of the 
solar convective zone provides time averaged constraints 
upon several physical parameters which enter the theoret- 
ical expression for the supply rate of energy, P, injected 
into the solar p modes. These constraints are: 

1) the depth dependency: of u 2 , the mean square 
velocity- of w 2 , the mean square vertical component of the 
velocity- of s 2 , the mean square values of entropy fluctu- 
ations. 

2) the wavenumber (fc) dependency of E and E s the 
turbulent kinetic energy spectrum and the turbulent en- 
tropy spectrum respectively. 

3) the depth dependency of the wavenumbers k^ and 
k^ s , the wavenumbers at which convective energy is max- 
imum and is injected into the turbulent inertial ranges of 
the turbulent kinetic energy spectra E, E a respectively. 



Differences between w 2 - and s 2 - and their respec- 
tive GMLT estimations have only small consequences on 
the profile of the excitation rate P{u)- However the val- 
ues reached by w 2 and s 2 with the 3D simulation are re- 
sponsible for an increase of the relative contribution of 
the Reynolds stress to P{oj) by a factor ~ 1.5 at low fre- 
quency v < 3 mHz compared to the one obtained with 
the GMLT solar model. This is because the GMLT model 
overestimates s 2 by ~ 20 % at the top of excitation region 
and underestimates w 2 within most part of the excitation 
region by up to ~ 15 %. 

The energy distributions E and E s over eddies with 
wavenumber k obtained in the simulation scale approx- 
imately as k +1 in the domain k < k^ . They therefore 
have approximately the same behavior as the 'Extended 
Kolmogorov Spectrum' (EKS) defined in Musielak et al. 
(1994). In contrast, their fc-dependencies significantly dif- 
fer from those assumed in the Nesis Kolmogorov Spectrum 
(NKS) which scales as k~ 5 below kg. The NKS predicts 
much larger maximum values for P than does the EKS. 
This is because the NKS concentrates kinetic energy in 
the vicinity of . 

The 3D simulation indicates that k$ ~ 3.6 Mm -1 at 
the top of excitation region. This corresponds to the hor- 
izontal size of the granulation (~ 2 Mm). It is worth not- 
ing that taking the depth dependency of ko(z) into ac- 
count results in a increase of the maximum of P by as 
much as ~ 50% and brings these values even closer to the 
observations. On the other hand, only minor differences 
are seen on the frequency dependence of the excitation 
rates P when using the depth dependency of fco( z ) from 
the simulation or assuming the form fc^ LT = 2tt//3A with 
A = aH p provided that j3 is adjusted in order for /c^ LT 
to match the value reached by kg at the layer where w 
reaches its maximum. 

The excitation rate P(to) is very sensitive to the value 
of $. In the GMLT formulation, the quantity $ is a pa- 
rameter which is adjusted in order to obtain the best 
fit between computed solar damping rates and the so- 
lar measurements: the adopted value is <& = 1.37 (see 
iHoudek et alJEoOlh . On the other hand, the 3D simula- 
tion suggests a higher value within the excitation region 
($ ~ 2). Larger values of $ result in an increase of the 
mode driving by the turbulent motions. We find that using 
the value $ = 1.37 underestimates P(oj) by a factor ~ 5 
relatively to P(u>) computed with &(z) in the excitation 
region from the 3D simulation. On the other hand, using 
the GMLT formulation for the convective velocity with a 
value $ ~ 2, as suggested by the 3D simulation, yields a 
power P(oj) close to the one obtained by the 3D simula- 
tion. To fix ideas, the maximum amplitude is ~ 4 cm/s, 
~ 8 cm/s, ~ 10 cm/s when calculated with GMLT and 
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$ = 1.37, with GMLT and $ = 2 and when using veloci- 
ties and $(z) derived from the 3 D simulation respectively. 
These figures must be compared to the observed maximum 
amplitude ~ 23 cm/s. 

This shows that the values of $ found for the so- 
lar GMLT model when adjusted to the damping rates 
is not compatible with the actual properties of the tur- 
bulent medium in the excitation region. An improvement 
could come from a consistent calculation which would as- 
sume a depth dependent $(z), as suggested by the simula- 
tions, in both damping rate and excitation rate computa- 
tions. Damping rates are indeed expected to be sensitive 
to depths deeper than the excitation rate where smaller 
values of $ are encountered and the simulation shows that 
the velocity anisotropy factor <E> decreases from 2 down to 
1.3 from top of the superadiabatic region to bottom of the 
simulated solar region. 

Without any adjustment of scaling parameters but us- 
ing all the constraints inferred from the 3D simulation 
considered here, we find a maximum of P much larger 
(~ 5 times larger) than the P maximum obtained using 
a ID GMLT solar model when <E> is fixed by the observed 
damping rates. It is also found that the so-called entropy 
source term, which arises from the advection of the tur- 
bulent fluctuations of entropy by the turbulent motions, 
is still the dominant source of the excitation. However its 
contribution to the excitation is now ~ 65 — 75% instead 
of ~ 95 % as found in Paper II. 

Our computation still underestimates by a factor ~ 3 
the maximum value of P compared with the one de- 
rived f rom the solar seismic observations by Chaplin et alJ 
l)l998(l . Moreover the decrease of P with v at high fre- 
quency (y > 3.5 mHz) is found to be significantly smaller 
than the one inferred from the solar seismic observations, 
indicating a deficiency in the present modelling at high 
frequency. 

As a final point, we discuss the model for the turbulent 
kinetic energy spectrum: 

In Paper II the parameter A and were adjusted 
- given a turbulent spectrum E(k) - so as to obtain the 
best possible agreement between computed and measured 
values of the maximum solar oscillation amplitude and its 
frequency position, as well as the frequency-dependence 
of the oscillation amplitudes. Adjustments of these scaling 
parameters led to a better agreement between computed 
values of P and the seismic observations when using the 
NKS than the EKS. However, the present results from a 
3D simulation strongly suggest that the EKS is a better 
model for the solar turbulent kinetic energy spectrum. 

The better agreement obtained with NKS than EKS 
when adjusting the free parameters is due to the fact that 
the NKS concentrates most of the kinetic energy in the 
vicinity of k^. Indeed, the NKS predominantly excites 
the modes whose period are close to the characteristic 
lifetime of the eddies of wavenumber fcjf, i.e. the modes 



with frequency close to the frequency at which P peaks 
(y ~ 3 mHz). As a consequence, the amount of energy 
going into the high frequency modes is relatively smaller 
with the NKS than it is with the EKS. This explains why 
the NKS reproduces better the steep decrease with v of P 
at high frequency and results in a value for k® identical 
to that inferred from the simulation (fcjf ~ 3.6 Mm -1 ). In 
contrast, whatever the adjustment, the EKS reproduces 
neither the frequency dependence of P at high frequency 
nor the value fcjf ~ 3.6 Mm" 1 . Hence assuming that the 
3D simulation yields the proper behavior of the solar ki- 
netic energy spectrum, well modelled by the EKS, one is 
led to conclude that the excitation as given by the present 
stochastic excitation model is not efficient enough at large 
scales (k ~ fco) and too efficient at small scales (k > fco). 

Discrepancies between our calcula tions and the ob- 
served excication rates or the results bv lStein k. Nordlundl 
(2001) are likely due to dynamic properties of turbulence 
which are not properly taken into account in the excitation 
model. Indeed, the dynamic properties of turbulence are 
modeled by the function Xk- All current theoretical cal- 
culations of th e excitation rates assume a gaussian func- 
tion for y fr fe.g.lGoldreich k, Keelevl <)l977l GK hereafter), 
iBalmforthl jl992a|) L The gaussian model is likely to be at 
the origin of the current under- estimate of the rates at 
which solar p-mode s are excited (see forthcoming paper 
ISamadi et alJliooii) . This may also explain the fact that 
we find that the entropy source term is dominant over the 
Reyno lds stress contribution whereas IStein k. Nordlundl 
(2001) in their direct computations found the reverse. In a 
recent study based on a frequency analysis of the present 
simulation we investigate what model can correctly repro- 
duce model Xk in the frequency range where the acoustic 
energy injected in to the solar p-mode s is important (see 
forthcoming paper ISamadi et al"ll2003|) . 



In the manner of [Rosenthal et all 



constraints 

from 3D simulation can be imposed to the ID model. 
According to the authors, such constraints result in a bet- 
ter agreement between the observed frequencies of the so- 
lar p-modes and the eigenfrequencies of the computed adi- 
abatic oscillations. An improvement in the calculation of 
the excitation rates at solar-type oscillations could then 
also come from a more consistent calculation of the eigen- 
modes which would use such constrained ID model. 
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